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ABSTRACT 

We present a new hybrid code for large volume, high resolution simulations of cosmic reionization, 
which utilizes a N-body algorithm for dark matter, physically motivated prescriptions for baryons and 
star formation, and an adaptive ray tracing algorithm for radiative transfer of ionizing photons. Two 
test simulations each with 3 billion particles and 400 million rays in a 50 Mpc//i box have been run to 
give initial results. Halos are resolved down to virial temperatures of 10 4 K for the redshift range of 
interest in order to robustly model star formation and clumping factors. This is essential to correctly 
account for ionization and recombination processes. We find that the halos and sources are strongly 
biased with respect to the underlying dark matter, re-enforcing the requirement of large simulation 
boxes to minimize cosmic variance and to obtain a qualitatively correct picture of reionization. We 
model the stellar initial mass function (IMF), by following the spatially dependent gas metallicity 
evolution, and distinguish between the first generation, Population III (PopIII) stars and the second 
generation, Population II (PopII) stars. The PopIII stars with a top-heavy IMF produce an order of 
magnitude more ionizing photons at high redshifts z > 10, resulting in a more extended reionization. 
In our simulations, complete overlap of HII regions occurs at z ~ 6.5 and the computed mass and 
volume weighted residual HI fractions at 5 < z < 6.5 are both in good agreement with high redshift 
quasar absorption measurements from SDSS. The values for the Thomson optical depth are consistent 
within 1 — a of the current best-fit value from third-year WMAP. 

Subject headings: cosmology: theory - large-scale structure of universe - intergalactic medium 
galaxies: formation - stars: formation - methods: numerical - radiative transfer 



1. INTRODUCTION 

Cosmic reionization starts when the first generation 
of stars begin to photoionize beyond the interstellar 
medium (ISM) and into the intergalactic medium (IGM). 
The overall picture is complex, with constant creation 
and possible demise, as well as mergers of individual 
HII regions. The mean ionized fraction will grow as the 
cosmic star formation rate increases. When a complete 
overlap of HII regions occurs, reionization comes to an 
end and the universe beco mes transparent to UV pho- 
tons. For r ecent r eviews, see lBarkana fc L ocb (200lf) and 
iFan et all (|2006af >. 

Currently, there are two major observational con- 
straints on the reionization epoch. SDSS observations 
indicate that reionization ends at z ~ 6 because spec- 
tra of quasars show lack of complete Gunn-Peterson 
absorption at lower redshifts, while complete Lyman 
alpha abs orption is found at immediately higher red- 
shifts (e.g. IFan et alJl25oH iBecker et ai]l2001t IFan et all 
2006b]). WMAP polarization measurements have yielded 
a Thomson optical depth of r = 0.09 ± 0.03, which sug- 
;ests that the IGM was largely ionized by redshift z ~ 10 
Page et al.ll20M ISpergel et al.ll2007h . The combination 
of these two observations imply that reionization may be 
considerably more complex than the simple assumption 
of an impulsive phase transition. Current and future ob- 
servations, primarily high redshift surveys of Lym an al- 
pha emitting galaxies fe.g. lKashikawa et al.ll2006l ). CMB 
measurements of the kinetic Sunyaev-Zeldovich (KSZ) 



effect from free electrons (e.g. ACT, SPT, Planck), and 
radio observations of the 21 cm radiation from neutral 
hydrogen (e.g. LOFAR, MWA, SKA), will provide sig- 
nificantly more detailed information. 

Progress on the theoretical front have mostly come 
from analytical and s emi-analytic al modelling (e.g. 
Barkana fc Loebl 12001 : JCenl l2003al lbl: IWvithe fc Loeb 
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simplified approaches, are able to afford continuous 
dynamic range and economically provide insightful in- 
formation. However, they are generally limited in scope 
and applicability. The semi-analytical approach provides 
relatively more realistic treatment of the formation of 

structure and HII region s. In the most recent work 

to da te (e.g. IZahn et all l2007t iMesinger fc Furlanettol 
|2007| ). realistic halo distributions and ionization maps 
have been generated by combining the excursion-set 
formalism with Lagrangian perturbation theory. These 
algorithms are computationally fast and are very useful 
for exploring the parameter space, especially when 
the prescriptions are calibrated in detail with realistic 
simulations. 

Numerical simulations directly solve the nonlinear 
physics of gravitational collapse, hydrodynamics, and 
radiative transfer. Furthermore, simulations allow a 
more straightforward and robust implementation of in- 
tricate prescriptions for astrophysical processes, such as 
star formation, which currently are not simulated from 
first principles. However, simulations of reionization are 
costly and hence, restricted in dynamic range because of 
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limitations in numerical algorithms and computational 
resources. High resolu t ion but small volume simula- 
tions (e.g. lGnediiJl2000t iRazoumov et al.ll2002ft can re- 
solve dwarf galaxies, Lyman limit systems, and mini- 
halos, allowing for more accurate determination of the 
star formation rates, photoionization rates, and clump- 
ing factors. However, they give highly biased results due 
to cosmic varia nce. Large volume but low resolution 
simulations (e.g. ISokasian et al.ll2003t iCiardi et al1 | 2003 |: 
Kohler et al.ll2005HHiev et al. l l2006aUbT]Zahn et alj|2007t 
McQuinn et al.ll2007f ) enablethe study of the large-scale 
structure of the dark matter, baryons, sources, and HII 
regions, but they must correct for relevant processes hap- 
pening on unresolved mass, spatial, and temporal scales. 

Until now, no other work reported in the literature has 
simultaneously satisfied the following two demands. On 
one hand, high resolution is required to resolve high red- 
shift halos with masses ~ 10 75 — 1Q s ' 5 Mq/1i, where the 
majority of photoionizing sources are thought to reside. 
Probing small scales is also essential to correctly calculate 
clumping factors and recombination rates. On the other 
hand, it has been argued that a large simulation volume 
of ~ (100 Mpc/h) 3 is necessary in order to fa irly sample 
highly biased sources (|Barkana fc Loebll2004h and large 
HII regions with chara cteristic sizes ~ 10 — 50 Mpc/h 
(|Furlanetto et al J 12004). The demand of having at least 
10 particles to identify a halo in a simulation with a vol- 
ume of (100 Mpc/h) 3 translates to a requirement of more 
than 20 billion particles. In this paper, we describe the 
methodology developed to cross this threshold of mass 
dynamic range and will present simulati ons with up to 24 
billion particles in companion papers (|Shin et al.ll2007l . 
Trac et al. 2007). 

Three-dimensional radiative transfer (RT) calculations 
in cosmological simulations are very challenging because 
the radiation field can be quite complex. By ray tracing, 
one can accurately propagate the radiation field, but it is 
computationally very expensive. Since the total number 
of rays emitted by a source scales with the size of the 
simulation N and the number of sources also scales with 
N, too many rays are generated with this 0(N 2 ) scaling. 
The prefactor is also relatively large because the number 
of sources can be as large as ~ 10% of N. We have 
developed an accurate and efficient adaptive ray trac- 
ing algorithm where w e use the ray splitting scheme of 
lAbel fc Wandel t (2002) to obtain excellent angular res- 
olution for ray tracing from point sources, but have in- 
troduced a new ray merging scheme in order to avoid 
the 0(N 2 ) scaling problem. With adaptive merging of 
co-parallel rays, the algorithm converges to O(N) scaling 
as the radiation filling factor approaches unity. 

In this paper, we use the cosmological parameters: 
tt m = 0.26, n A = 0.74, n b = 0.044, h = 0.72, ct 8 = 0.77, 
and n s — 0.95, based on the lat est results from WM AP, 
SDSS, BAO, SN, and HST (see lSpereel elaTl [20071 and 
references therein). $2] describes the n umerical methods : 
a N-body particle-multi-mesh (PMM; iTrac fc Per3l2006D 
code for the gravitational evolution of the dark matter, 
a hydrostatic equilibrium model for gas in halos, a star 
formation prescription for gas cooling in halos, and a 
RT algorithm with adaptive ray tracing. fJJ] presents two 
high-resolution simulations each with 3 billion particles, 
400 million rays, and 180 3 RT grid cells in 50 Mpc/h 
boxes. In one simulation, we consider only Population II 



(PopII) stars with a Salpeter initial mass function (IMF), 
while in the second simulation, we include Population III 
(PopIII) stars with a top-heavy IMF. if4] discusses initial 
results: halo mass functions, halo clustering, star for- 
mation rates, ionization fractions, clumping factors, and 
optical depths. We demonstrate that the current ob- 
servational constraints require a star formation history 
where both PopIII and PopII stars are significant sources 
of photoionizing radiation. 

2. NUMERICAL METHODS 

2.1. Dark matter 

We use a particle-multi-mesh (PMM; ITrac fe Pen! 
2006() N-body code, based on an improved particle-mesh 
(PM) algorithm, to compute the gravitational dynam- 
ics of the collisionless dark matter. In principle, stan- 
dard PM codes can achieve high spatial resolution with 
a large mesh but this comes at a great cost in memory 
and to a lesser extent in work. In practice, they are nor- 
mally limited to a mean interparticle spacing to mesh cell 
spacing ratio of 2:1, where the storage requirements for 
particles and grids are approximately balanced. PMM 
utilizes a domain-decomposed, FFT-based gravity solver 
to achieve higher spatial resolution without increasing 
memory costs. It is based on a two-level mesh Poisson 
solver where the gravitational forces are separated into 
long-range and short-range components. The long-range 
force is computed on the root-level, global mesh, much 
like in a PM code. To achieve higher spatial resolution, 
the domain is decomposed into cubical regions and the 
short-range force is computed on a refinement-level, local 
mesh. In the current version, PMM can achieve a spatial 
resolution of 4 times better than a standard PM code at 
the same cost in memory. The simulations in this paper 
were run with a mean interparticle spacing to mesh cell 
spacing ratio of 4:1. 

For each N-body particle of fixed mass, we store its 
comoving position x and velocity v in order to solve 
the equations of motion in an expanding Friedman- 
Robertson- Walker (FRW) universe. In addition, the lo- 
cal matter density p m associated with each particle is 
calculated using a SPH-like smoothing kernel. These 
seven particle variables are updated every PM time step 
of AtpM < 10 7 years. 

We locate dark matter halos using a friends-of-friends 
algorithm with a standard linking length of 0.2 times the 
mean interparticle spacing. Halo catalogs are produced 
containing the following information: redshift, comoving 
positions and velocities, and mass. In addition, we tag 
each particle that is bound to the identified halos in order 
to model baryons and star formation. 

In a SCDM cosmology with fl g = 1 and Vi\ = 0, the 
spherical top-hat collapse model ([Pee bles 1980]) defines a 
virialized halo to be a sphere of mass M v - lr and radius i? V i r 
with a characteristic average densit y equal to A e = 18-7T 2 
times the critical density p cr it(z) (|Lacev fc Cold [l993l 
119941) . In a ACDM cosmology with tt A ^ 0, the char- 
acteristic average density for virializa tion has been cal- 
culat ed using numerical simulations (|Brvan fc Normanl 
1998), but note that at the high rcdshifts z > 6 prior to 
complete reionization, the universe is basically Einstein- 
de-Sitter and A c ss 187r 2 . 

2.2. Baryons 
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On large, linear scales, we assume the baryons are unbi- 
ased tracers of the dark matter because the hydrodynam- 
ics of the baryons is primarily dictated by the large-scale 
gravitational dynamics of the matter and not by the ther- 
modynamics of the gas. On small, nonlinear scales, the 
collisionless dark matter continues to collapse, undergo- 
ing violent relaxation and forming virialized halos with 
deep potential wells. The collisional baryons fall into the 
potential wells, get shock heated to approximately the 
virial temperature, and the collapse gets halted by gas 
pressure. However, this collisional component can dis- 
sipate energy through radiative cooling and eventually 
collapse to very high densities to form stars. 

We approximate the local baryon density pt, associated 
with each particle by multiplying the matter density p rn 
by the cosmic baryon fraction Q;,/f2 m , except in high 
density regions where gravitational collapse and pressure 
support are important. In virialized halos, we calculate 
the gas density profile assuming the gas is in hydrostatic 
equilibrium with the gravitational potential. 

For a virialized halo with mass M v ; r and radius i? v ir, 
the dark matt er density profile is taken to follow the 
NFW formula (jNavarro et al.lll997| ): 



Pdm(x) = 



M vir c 3 /(4^. r ) 



1 



ln(l + c) - c/(l + c) x(l + x) 2 



(1) 



where x = r/r s , r s is a scale radius, and c = R v - n -/r s is a 
concentration parameter. Simulations of early structure 
formation in a ACDM cosmology produce dark matter 
halos which are well-fitted by the NFW formula for red- 
shifts z < 12 a nd with only sma ll deviations for redshifts 
up to z ~ 49 (|Gao et al.ir2005l ). The concentration pa- 
rameter is approximately given by (see lDolag et~aTll2004l . 
and references therein), 



c(M vil ,z) 



7.5 



1 + z \ M, 



(2) 



where the nonlinear mass is M* = 2.09 x 10 12 M Q //i for 
the chosen cosmology. The high resolution of the sim- 
ulations is sufficient to locate halos and measure their 
masses, but it is not enough to accurately measure their 
dark matter profiles. Therefore, the NFW profile and 
the parametric model for the concentration parameter 
are used instead. 

Assuming hydrostatic equilibrium and a polytropic 
equation of state P g oc pj, the gas density and tempera- 
ture profiles can be parametrized as, 



Pg( x ) = PoVg(x) , 
T g {x)=T y1-\x) 



(3) 
(4) 



where 7 is the polytropic index and the coefficients 
po = p g (0) and To = T g (0) are two boundary conditions. 
The dimensi onless gas density profile has the analytical 
solution (e.g. [Komatsu & Seljak 200j]), 



y g (x) = il-A 



1 - 



ln(l + x) 



1/(7-1) 



where 



A = 



1\ ( GM vil n 



ln(l 



(5) 



(6) 




Fig. 1. — The mass range of halos where star formation can take 
place. The curve M(z, /) gives the mass range M > M(z, /) which 
accounts for the fraction / of the total star formation at redshift 
z. The five solid curves from top to bottom are for the fractions 
/ = (0.2, 0.4, 0.6, 0.8, 1.0). For comparison, the simulations all 
have a particle mass of m p = 3.02 X 10 6 Mq/A (dotted line) and a 
minimum halo mass of Mmin — 20m p , (dashed line). 

and p is the mean molecular weight. 

We have chosen 7 = 1.2 as suggested by numerical 
simul ations and semi-analytic models (see lOstriker et all 
I2005L and references therein). The baryon fraction at 
the virial radius is normalized to be equal to the cosmic 
value. We also normalize the temperature at the virial 
radius to be equal to the halo virial temperature, 



T ■ — 



GM V 



3kR v 



3k 



1/3 



[GM vir H(z)] 



2/3 



(7) 



Given a dark matter density, we first analytically solve 
a cubic equation to find the radius and then the corre- 
sponding baryon density. This is done for every parti- 
cle identified as belonging to a halo and for every PM 
timestep. 

2.3. Star formation 

Halos with virial temperatures > 10 4 K are thought 
to form stars more efficiently than mini-halos. Above 
this transition temperature, efficient atomic line cooling 
allows the gas to dissipate energy and collapse to high 
enough densities for molecular hydrogen (H2) formation. 
Molecular cooling then allows the gas to further collapse 
to the very high densities needed for star formation in 
giant molecular clouds. Star formation in mini-halos is 
possible because of H2 formation and molecular cooling 
in dense self-shielded gas. However, the H2 in mini-halos 
are very susceptible to dis sociation by Lyman - Werner 
photons from the first stars (|Haiman eHdll997h . There- 
fore, we will assume that mini-halos make only minor 
contributions to the total stellar density and we model 
star formation only in halos which cool through atomic 
transitions. 

In our simulations, the star formation criteria re- 
semble those used in hydrodynamic simulations (e.g. 
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ICen fc Ostrikeijri992f h The first criterion for star forma- 
tion is that it only occurs in halos with virial tempera- 
tures above 10 4 K. The second criterion is that the local 
cooling time must be shorter than the local dynamical 
time: 

3p g kT 



^cool — 



^dyn 



2// 



K(A-r) 



3tt 



32Gp„ 



(8) 



(9) 



(A-r 



where X is the hydrogen fraction by mass and n 2 M L 
is the net cooling rate per unit volume (see iCenl 1199 
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iTheuns et al.l H998). We calculate the cooling radius as 
a function of virial mass and redshift assuming that the 
gas in the halos are in ionization equilibrium, which is a 
good approximation because of the high densities. Feed- 
back from photoionization and supernova will affect the 
cooling radius, but modeling these affects is beyond the 
scope of this paper and will be left for future work. Fig- 
ure [1] shows the mass range of halos where star formation 
can take place. The curve M(z, /) gives the mass range 
M > Al(z,f) which accounts for the fraction / of the 
total star formation at redshift z. The five curves from 
top to bottom are for the fraction / = (0.2, 0.4, 0.6, 0.8, 
1.0) of the total star formation rate. 

For every particle within the cooling radius of a T V1I > 
10 4 K halo, we calculate its star formation rate as 



dt 



C*Mg 

id 



(10) 



where M g = Mb — Af* is the gas mass and c* is a star 
formation efficiency. We keep track of the stellar frac- 
tion /* for each particle in order to differentiate the 
total baryonic mass into gas and stars. Furthermore, 
we distinguish between the first generation, PopIII stars 
and the second generation, PopII stars. It is postu- 
lated and this is supported by hydrodyna mic simulations 
(jO'Shea fc Norman! |2005 lYoshida et al] I2006D that the 
first generation of stars forming from pristine primordial 
gas are very massive M ~ 100M© and the stellar IMF at 
high redshifts is top-heavy compared to a Salpeter IMF. 

We have a simple prescription for the formation of 
PopIII stars motivated by results from hydrodynamic 
simulations. Once a halo has accreted enough mass such 
that its virial temperature exceeds 10 4 K, star formation 
is turned on and the IMF is determined by the mctallicity 
of the halo. An initially metal-poor halo with metallicity 
less than some critical metallicity Zp op iii first undergoes 
PopIII star formation for a duration of ip op ni years af- 
ter which it switches to PopII star formation. Hydrody- 
namic simulations have found that once the metallicity 
reaches a value ^p op iii = 10~ 35 Zq, metal line cooling 
becomes efficient enough that giant molecular clouds can 
undergo more fragmentation and form less massive stars. 
The halo self-enrichment time ip op in = 20 Myr takes into 
account the lifetime (~ 3 Myr) of a massive star and the 
time it takes for supernova feedback to traverse and en- 
rich the interstellar medium (ISM) with metals produced 
by the stars. Furthermore, a fraction /z, C sc = 10% of 
the metals are assumed to escape into the intergalactic 
medium (IGM) and we propagate them assuming a con- 
stant average velocity i>z,igm = 10 km/s. By tracking 
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Radiation filling factor 

In the adaptive ray merging scheme, 



maximum number of rays with levels I > i mer ge in a cell, but the 
total number of rays can be larger due to contributions from rays 
at lower levels, which remain singular. In the limit of complete 
reionization, the number of active rays will converge to a value of 
approximately N meTge times the number of RT grid cells and the 
RT algorithm will scale as O(N). 

the metals, we can suppress PopIII star formation us- 
ing a local metallicity rather than relying on a global, 
volume-averaged metallicity for the IGM. 

For PopII stars with a Salpeter IMF, the number of 
ionizing photons per baryon of star formation is 5200, 
4100, and 270 for the three frequency ranges v > 13.61, 
24.59, and 54.42 eV, respectively. For PopIII stars with 
a top-heavy IMF, the correspond i ng nu mbers are 70000, 
55000, and 3500 (jSchaereij 12001 |2003| ). The radiation 
escape fraction / 7jC sc is defined to be the fraction of ion- 
izing photons which escapes the ISM. 

2.4. Radiative transfer 

Our unique approach to three-dimensional radiative 
transfer (RT) complements the high resolution and ac- 
curacy of the N-body simulations. We have developed 
an adaptive ray tracing algorithm to accurately and effi- 
ciently propagate the complex radiation field generated 
by the many sources. Furthermore, we calculate ion- 
ization and recombination for each individual particle, 
taking into account their geometrical cross-sections and 
volumes in order to correctly account for the clumping 
factors and self-shielding. In all previous large volume 
simulations of reionization, the RT is calculated using 
low resolution density fields defined on a coarse grid, but 
small-scale information is lost with this approach. 

First, we describe the set up for the RT calculations. 
Particles, sources, and rays are collected on a grid with 
N-rt cells, where the mean number of particles per RT 
cell is 8 3 and the number of PM cells per RT cell is 32 3 . A 
source cell in the RT grid has a total star formation rate 
equal to the sum from all star-forming particles within 
the cell. A radiation cell has a total number of photons 
equal to the sum from all rays intersecting it. Each ray 
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carries with it three numbers of photons for the frequency 
ranges (eV): 13.61 < v < 24.59, 24.59 < v < 54.42, and 
54.42 < v < oo, respectively. Particles are treated as 
individuals, but their differing photoionization rates will 
depend on the photon number density of the cell in which 
they belong to. The HI, Hel, and Hell fractions for each 
particle are kept tracked of. The radiative transfer time 
step AtRT is set by the light-crossing time for a RT cell 
and there are ~ 10 — 100 RT time steps per PM time 
step. 

We have implemented an adaptive ray t racing algo- 
rithm similar to that of lAbel fc Wandeltl (|2002f) . but 
have introduced a new merging scheme in order to avoid 
the 0(N 2 ) scaling problem. The ray tracing algorithm 
makes use of the HEALPix (|G6rski et all 120021 l2005f) 
equal-area pixelization scheme, which ensures that the 
rays uniformly cover the unit sphere. We start with a 
base level Zbase = 1 and cast Abase = 12 x 4 ibaac = 48 
rays per radiation source cell each radiative transfer time 
step. In order to avoid artifacts arising from discreteness 
effects, one can randomly rotate the local coordinate sys- 
tem centered on each source cell, but we find it simpler 
and just as effective to randomize the origin of the lo- 
cal coordinate system within each source cell every time 
step. Once a source cell becomes transparent, its 26 near- 
est neighbours will each, on average, have approximately 
1.8 rays with radius r = 1 RT grid unit. We find that us- 
ing a higher base level Zbasc = 2 and casting Abase = 192 
rays is unnecessary because of the randomization within 
source cells, the large number of RT times step per PM 
time step, and the large number of sources contributing 
to the radiation field. 

As rays propagate, they adaptively split to ensure that 
a minimum number of rays, iV m i n , from the same source 
intersect each radiation grid cell in the path each time 



where N„ 



is the maximum number of rays and iV rac j 



is the number of radiation cells intersected by at least 
one ray. The ratio N Tayjmax /NnT sets the value we want 
in the limit of complete reionization and then it is di- 
vided by the radiation filling factor. Eq. (TT2"|) ensures that 
the number of active rays iV ray will converge to Af ray!inax 
rather than growing quadratically with the size of the 
simulation. 

The adaptive ray merging scheme is also implemented 
within the framework of HEALPix. For each RT time 
step, we define the level for merging as, 



merge 



Integer 



log(AT ray/RT /12) 



log(4) 



(13) 



where the Integer operator can round down, to the near- 
est, or up. Rounding down is the conservative approach 
to satisfying Eq. (|12|). but it unnecessarily lowers the an- 
gular resolution when the radiation filling factor is small. 
Rounding up increases the angular resolution, but it may 
violate Eq. (TT2"]) when the radiation filling factor is close 
to unity. A good compromise, while still maintaining sta- 
bility, is to round to the nearest integer. For each cell, 
only rays with levels / > l meT ge are considered for merg- 
ing, while ones at lower levels remain singular. The unit 



sphere is subdivided into N n 



12 x 4 



equal- 



area pixels and rays with I > Z me rge are subdivided into 
-Emerge bins. Rays with bin number Zb m are considered 
approximately parallel if they originate from rays of level 



I = I 



merge 



which share the same pixel number i plx = ibii 



Rays within a bin are merged into one and the total 
number of photons is conserved. The new position, di- 
rection, and effective radius are calculated by taking av- 
erages weighted by the photon count of the individual 
rays. 

In Figure the number Af mer gc is plotted as a func- 



step. A ray with level I and pixel number i pix splits tion of the radiation filling factor for N Tay / aT = 64. Note 



into 4 daughter rays with levels I 



4i p ix + n where n runs from to 3 



1 and pixel numbers 
We use a value of 

Nmin = 1-8 to be consistent with the base numbers for 
sources and this value is more than sufficient by the very 
same arguments given above. A ray of radius r in RT 
grid units will have a level, 

log(AU„^ 2 /3) 



l(r) = Integer 



(11) 



log(4) 

and since the Integer operator is always taken to round 
up, most cells will actually be intersected by more than 
N m in rays from the same source each time step. The 
adaptive splitting scheme provides excellent angular res- 
olution for ray tracing from point sources, but it is com- 
putationally expensive because of the 0(N 2 ) problem. 

We have developed an adaptive ray merging scheme 
which converts the initially 0(N 2 ) scaling to 0(A) as 
the radiation filling factor approaches unity. In our sim- 
ulations, source cells can occupy up to ~ 10% of the RT 
grid and therefore many rays will intersect any given cell 
at any given time. Therefore, as the radiation filling fac- 
tor increases, the angular resolution of rays entering and 
exiting cells can decrease without sacrificing much accu- 
racy. In our ray merging scheme, the number of rays in 
a cell at any given time is adaptively restricted with the 
simple formula: 



N, 



ray/RT 



ray .max 



A, 



RT 



A, 



RT 



rad 



(12) 



that A mergo is the maximum number of rays with levels 
I > Emerge in a cell, but the total number of rays can be 
larger due to contributions from rays at lower levels. For 
the majority of the simulation where the radiation filling 
factor is < 2/3, Amerge — 4iVbase and the angular resolu- 
tion in any given cell is better than the angular resolution 
of base rays in the source cells themselves. We have ex- 
perimented with higher angular resolution by increasing 
^base to 2 and N ray /in to 256, which comes with the ex- 
tra cost of a factor of 4 in both computational work and 
memory, but find no significant effects on the radiation 
field. 

Rays are collected on the grid each RT time step using 
the nearest grid point (NGP) assignment scheme. This 
mapping is robust because of the large number of rays in- 
tersecting any given cell and the large number of RT time 
steps per PM time step. Our ray tracing algorithm share s 
several similarities to that of McQuin n et al.l (1 2007). 
who modified a previous algorithm (jSokasian et al . 2001, 
2003). The algorithms have several commons, including 
how rays are cast, split, and mapped to the grid, but 
the main difference is how the problem of too many rays 
is solved. They consider a simple approach where rays 
no longer split after traveling a certain distance. The 
critical distance is chosen adaptively depending on the 
relative luminosity of the ray to the total luminosity of 
all rays in the box. However, this isotropic value is not 
optimal for HII regions, which are quite anisotropic. Our 



G 



Trac & Cen 



adaptive ray merging scheme is a more general approach 
and the ray tracing algorithm can be applied to many 
problems in astrophysics. 

Once the photon density field is updated on the grid, 
we calculate the number of ionizations and recombina- 
tions for each individual particle rather than using the 
cell-averaged hydrogen and helium densities. This ap- 
proach maximizes the resolution and prevents the clump- 
ing factors from being severely underestimated. Further- 
more, working with the particles directly gives us the ad- 
vantage of being able to account for self-shielding. Star- 
forming halos which are sources will be ionized from the 
inside out, but the structure which are sinks will gen- 
erally be ionized from the outside in. To allow for self- 
shielding, the particles within each RT cell are sorted 
from lowest to highest density and the photons are ab- 
sorbed by particles in this order. 

The equations of ionization evolution for hydrogen and 
helium are: 



dn 



HI 



dt 



dn 



Hel 



dt 



dnm 



in 



dt 



= affli"-c«Hi - Thi^hi, (14) 



= aHciin e nHeii - rnei^Hci, (15) 



'Q;HcIII«c^HcIII + ThcHI-HcII, (16) 



where a are the recombination rates and T are the pho- 
toionization rates. Consider a particle with volume V p , 
cross-section A p , and length l p in a RT cell where the 
number of photons per unit volume per unit frequency 
at frequency v is r\ v (cm -3 Hz -1 ). The rate of change of 
neutral hydrogen in the particle due to photoionization 
is given by 



^7Hi = Thi^-hi 



-dv, (17) 



where N v = r\ v A v c (s _1 Hz -1 ) is the number of pho- 
tons per unit time per unit frequency passing through 
the particle, tvhi = "hioVhi^p i s the optical depth for 
absorption by neutral hydrogen, and the weight, 



«VHI 



HHI0VHI 



HHICi/HI + "HclOVHcI + ^HcIIfi/HcII 



(18) 



takes into account the competition between HI, Hel, and 
Hell for the same photons. Source particles are set to be 
fully ionized since the radiation escape fraction is defined 
to be that wh ich escapes t he IS M. Eq. (| 1 T[) is similar to 
Eq. (15) from I Abel et al.l (|1999f ). except we have gener- 
alized it to account for a range in frequencies. The rates 
of change for Hel and Hell by photoionization similarly 
follow that for HI. 

Since the equations of ionization equilibrium are stiff, 
the time integration is computed using a stable implicit 
scheme rather than an explicit one. Time integration is 
performed over the RT time step, but often subcycling 
in steps of the recombination time is required since the 
latter time scale is generally smaller. The rays have their 
photon numbers reduced in proportion to the fraction of 
consumed photons throughout the cell. In a cell where 
< 0.01% of the photons remain unconsumed, we delete 
all the rays within in order to be computationally efficient 
and this has negligible effect on the results. 



3. SIMULATIONS 

The hybrid N-body plus RT simulations of reioniza- 
tion were run with the following cosmological parame- 
ters: n m = o.26, n A = o.74, n b = 0.044, h = 0.72, 

(78 = 0.77, and n s = 0.95. The N-body initial conditions 
were generated for an initial redshift z — 300 where the 
matter power spectrum is still linear and the real space 
density field has |<5 max | < 1 and as ~ 0.03. The initial 
matter transfer function wa s computed with CMBFAST 
(|Seliak fc Zaldarriagalll996T ) and the Zeldovich approxi- 
mation was used to calculate the displacement and ve- 
locity for the particles. 

Table Q] lists the parameters for the two test simula- 
tions, each with 3 billion particles in a 50 Mpc/h box. 
For the PM calculations, the comoving grid spacing is 
Axpm = 8.68 kpc/h and the particle mass resolution is 
m p = 3.02 x 10 6 Af Q //i, which is a substantial improve- 
ment over all previous large-scale simulations of reioniza- 
tion. Halos are identified down to a minimum mass of 
Mm ^ = 6 x 10 7 M Q //i. The RT grid has a comoving grid 
spacing of Accrt = 278 kpc/h and sources and photons 
are collected at this resolution. However, particles are 
treated individually for the ionization and recombination 
calculations, allowing us to probe the clumping factors 
down to scales of Axpm rather than Azrt- Figure [3] is 
a sample visualization of the HII density evolution with 
PopIII stars. 

The simulations were run at the National Center 
for Supercomputing Applications (NCSA) on a shared- 
memory SGI Altix with Itanium 2 processors. Each test 
simulation used 128 processors and 360 GB of memory 
and required approximately 12000 cpu hours to complete. 
Furthermore, we have completed a large simulation with 
24 billion particles in a 100 Mpc/ft, box. The LlOOa sim- 
ulation, run with an alternative RT algorithm, used 512 
processors, 2 TB of memory, and approximately 80000 
cpu hours. Detailed results will be presented in a com- 
panion paper (Shin et al. 2007). 

4. RESULTS 

4.1. Halo mass functions 

We identified dark matter halos using a friends-of- 
friends (FoF) algorithm with a standard linking length 
of b = 0.2 times the mean interparticle spacing. Fig- 
ure @] shows the differential mass function dn/dlog(M), 
where n(M, z) is the comoving number density of halos 
with mass less than M at redshift z. Figure [5] shows 
the cumulative mass function N(M,z), defined as the 
number of halos with mass greater than M at redshift 
z. The FoF res ults from the L50 sim ulations are in best 
agreement with IWarren et aTl (|2006l ). who conducted a 
careful measurement of the mass function of dark matter 
halo s with a suite of hi gh-res olutio n N-body simulations. 
The ISheth fc Tormenl (ST; 119990 model also provides 
a good fit, but still differs by several tens of percents, 
with larger differen ces at high er redshifts. The original 
iPress fc Schechterl (PS; I1974T ) model underpredicts the 
halo abundance at all our redshifts. Our results are con- 
sistent with recent work on the mass function of high 
redshift dark matter halos (jReed et al.ll2007t iLukic etafl 
[20071 : iCohn fc Whitdl2007h and arc further confirmation 
that standard PS underpredicts the abundance of high 
redshift halos, particularly at the high mass end. 
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TABLE 1 
Simulation parameters 



Model 


L (Mpc/ft) 


Np 


Npm 


Af RT 


Nray ,max 


Stars 


Comments 


L50a 


50 


1440 3 


5760 3 


180 3 


(1440 3 )/8 


PopII 


/ 7 ,GSC = 0.16 


L50b 


50 


1440 3 


5760 3 


180 3 


(1440 3 )/8 


PopII & PopIII 


/ 7 ,esc = 0.15, tp piii = 20 Myr, Uz.IGM = 10 km/s 


LlOOa 


100 


2880 3 


11520 3 


360 3 


N/A 


PopII & PopIII 


Shin ct al. (2007) with alternative RT algorithm 


LlOOb 


100 


2880 3 


11520 3 


360 3 


(2880 3 )/8 


PopII & PopIII 


Trac et al. (2007) 




Fig. 3. — Low resolution image showing the HII density in 1 Mpc/ft deep slices through the L50b simulation box at redshifts z = 9, 8, 
7, and 6. The volume-weighted mean HI fraction are /hi = 0.67, 0.46, 0.12, and 1.5 X 10 -4 , respectively. Higher resolution images can be 
found at http:/ /www. astro. princeton.edu/~htrac/reionization. html 
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Fig. 4. — Friends-of-fricnds differential mass functions of dark 
matter halos. The halos are identified using a standard linking 
length of 0.2 times the mean interpa r ticle s pacing. The simulation 
results are best fit by | Warren et al.l lj200S. blue, solid) at all red- 
shifts. ISheth &: To rmcn ( 1993, green , dashed) tends to give h igher 
abundances at higher redshifts, while Press & Sclicchtcr (1974, red, 
dotted) underpredicts the abundances at all redshifts. 
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Fig. 5. — Friends-of-fricnds cumulative mass functions of 
dark matte r halos at z = (6, 9, 12, 15) with Poisson error 
bars. The Warren ct al. (2006, blue, so l id) p rediction provides 
a better fit than both Press fc Schechterl 1 119741 . red, dotted) and 
ISheth fc Tormenl ^19991 . green, dashed). 

In earlier preliminary work, we used simulations of the 
same size but with a late starting redshift of z = 60. 
We also previously used a spherical overdensity (SO) al- 
gorithm with a characteristic density chosen to be 200 
times the crit i cal de nsity, similar to that described in 
lLacev fc Colel (|1994l) . The SO results appeared to be 
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Fig. 6. — Dark matter halo power spectra are shown for five log- 
arithmic mass (MQ/h) bins (from bottom to top): 8.0-8.5 (red), 
8.5-9.0 (orange), 9.0-9.5 (green), 9.5-10.0 (blue), and 10.0-11.0 
(magenta). The nonlinear (black, solid) and linearly extrapolated 
(black, dashed) matter power spectra are also shown for compari- 
son. 

in better agreement with PS, particularly at higher red- 
shifts and higher masses. However, we were systemati- 
cally under-resolving halos because our starting redshift 
was too late to a ccurately capture the nonlinear gravi- 
tational collapse. iReed et"aT1 (|2007| ) have suggested that 
simulations must start ~ 10 — 20 expansion factors before 
the redshift at which results converge. Another system- 
atic effect came from our SO definition of halo mass. 
M200 is generally smaller than MFoF and a mass func- 
tion constructed with the former will h ave a lower am- 
plitude than the latter at a fixed mass. iCohn fc White] 
(|2007f ) have quantified the difference in the mass function 
when two FoF linking lengths of 0.2 and 0.168 are used. 
While the mass function for the larger linking length is in 
good agreement with ST, the results for the smaller link- 
ing length are better fit by PS but still with discernible 
differences. 

4.2. Halo clustering 

Bar kana fc Loebl (|2004( ) have emphasized that the high 
redshift halos where photoionizing sources reside are 
highly biased and this is clearly seen in our simulations. 
In Figure [6l we plot the dimensionless halo power spec- 
trum, 

A 2 hh (k,M,z) = A^(|,5 h (fc,M,z)| 2 ), (19) 

for wavenumbers k where the signal to noise is greater 
than unity. The Poisson noise is subtracted by removing 
the white noise power due to all self-pairs. All halo spec- 
tra resemble power-laws, regardless of mass and redshift, 
and there appears to be an inverse relationship between 
the effective slope and the halo number density. This 
suggests that the nonlinear clustering of dark matter ha- 
los can be described with a self-similar parametrization. 
We will quantify this relationship in upcoming work. 
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Fig. 7. — Dark matter halo bias are shown for five logarithmic 
mass (Mq/K) bins (from bottom to top): 8.0-8.5 (red), 8.5-9.0 
(orange), 9.0-9.5 (green), 9.5-10.0 (blue), and 10.0-11.0 (magenta). 
The bias increases with mass and with redshift and is only linear 
on scales where.^ is small. 

In Figure u\ we plot the halo bias, 



b(k,M, z) = 



Ay 

111 



(k,z) 



(20) 



where A^ in (A;, z) is the linear matter power spectrum ex- 
trapolated to redshift z. The halo bias increases with 
mass and with redshift and is only linear on scales where 
#h is small. The linear bias h as been derived for the 
PS model (IMo fc White! [l99^ and for the ST model 
(|Sheth fc Tormenl 119991 7 PS predicts higher bias than 
the simulations because in this model, halos are more 
rare at a fixed mass and redshift. On the contrary, the 
ST prediction gives slightly lower values than our results. 

Since the linear bias is inversely related to the halo 
number density, the derivative dfr/dM increases rapidly 
with mass because of the exponential decline in the abun- 
dance of massive halos. This rapid change is particu- 
larly more pronounced at higher redshifts. The mass- 
dependent bias of halos has important impli cations for 
the clu stering of sources and HII regions. Illiev et ail 
(2006aj) were only able to resolve halos with M > 2.5 x 
10 9 M Q , which is a factor of 30 times larger than our min- 
imum halo mass. In their simulation, reionization ended 
at z w 12 and at this redshift, their power spectrum of 
all resolved halos is approximately a factor of 3 larger 
than ours. The relative bias will be even larger at higher 
redsh ifts. Recent work (Zahn et al. 2007; McQuin n et al.l 
|2007|) have pointed out that this leads to a very different 
picture of reionization, characterized by relatively fewer 
but larger and more spherical HII regions. Therefore, 
simulations of reionization must resolve halos down to 
a minimum mass where the virial temperature is ~ 10 4 
K in order to correctly capture the clustering of sources 
and HII regions. 

4.3. Star formation 




Fig. 8. — The top panel shows the comoving star formation 
rate (SFR) from the L50b simulation with PopIII stars. The total 
SFR (black, solid), from PopIII (red, long-dashed) and PopII (blue, 
short-da shed) stars, is in good agreem ent with the semi-analytical 
model of Hernquist & Springcl (2003, black, dotted). The bottom 
panel shows the cumulative stellar density as a fraction of the mean 
baryon density. 




Fig. 9. — The top panel shows the comoving photon produc- 
tion rate d(n 7 /riH)/dt of photons with energies > 13.61 eV. The 
total rates from the L50a and L50b simulations are given by the 
lower and upper (black, solid) curves, respectively. In the L50b 
simulation, the PopIII (red, long-dashed) stars contribute a larger 
fraction to the total rate than the PopII (blue, short-dashed) stars 
for z > 12. The bottom panel shows the cumulative photon density 
as a fraction of the mean hydrogen density. In the L50b simulation, 
the photon budget is dominated by PopIII stars for z > 10. 
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In principle, the cosmic star formation rate (SFR) 
can be modelled using hydrodynamic simulations 
(e.g. ISpringel fc Hernquistl [2001 iNagamine et alj 12004ft 
and calibrated aga i nst l ow redshift observations. 
iHernquist fc Spring cl (2003) have derived an analytical 
formula for the SFR by considering the mass within the 
cooling radius of collapsed halos. The free parameters 
of the model are fitted using results from hydrodynamic 
simulations. In our simulations, the shape of the curve 

(z) agrees well with theirs and we choose the star for- 
mation efficiency c* = 0.03 to match their normalization, 
corrected for our cosmology. It is encouraging that our 
star formation history is very similar to that obtained 
in hydrodynamic simulations and semi-analytic models. 
Note that in our simulations, the number of escaped ion- 
izing photons depends on the product c*/ 7iesc . While the 
star formation efficiency and the radiation escape frac- 
tion are degenerate in this regard, it is still important to 
model their values individually since other calculations 
will require different combinations of these two param- 
eters. For instance, in scatt er calculations fo r high red- 
shift Lyman alpha emitters (|Tasitsio mi 2006), the prod- 
uct c*(l — / 7 , C sc) is also important. 

Figure [5] shows the comoving SFR and the cumulative 
stellar density as a function of redshift. The L50a and 
L50b simulations have the same total SFR, but the lat- 
ter differentiates between PopII and PopIII stars. While 
the stellar density of PopII stars is generally larger than 
that of PopIII stars, the number density of photoionizing 
photons produced by the former does not exceed that of 
the latter until z — 10, as shown in Figure [5J We assume 
that PopIII stars have a top-heavy IMF and produce a 
factor of 13.5 more photoionizing photons per unit stel- 
lar mass. In the L50b simulation, the photon produc- 
tion rate dn 7 /di changes more gradually with redshift. 
We show that this leads to a more extended reionization 
epoch in §4.41 and to a larger Thomson optical depth in 
H4.6I when the redshift of complete reionization is fixed. 

The SFRpopni reaches a broad maximum value of 
~ 10~ 3 M Q /Mpc 3 /yr at z ps 1 0, which is consis- 
tent with semi-anal ytical models (jYoshida et alJ 120041 : 
IWvithe fc Cer3l20"06| ). We find that the decline of PopIII 
stars at lower redshifts is mainly due to the following rea- 
son. In our simulations, PopIII stars form in metal-poor 
halos which have just accreted enough mass such that 
the virial temperature exceeds 10 4 K. However, the cool- 
ing mass M v ; r (T v i r = 10 4 K) is increasing as the redshift 
decreases (see Fig. [J) and the formation rate of metal- 
poor halos with this mass is decreasing in time. We also 
find that the amplitude of the SFRp op ni is more strongly 
dependent on the halo self-enrichment time tp op iii than 
on the metal enrichment of the IGM, which is controlled 
by the escape fraction of metals /z, C sc and the average 
propagation speed of metals vzjgm- The relative contri- 
bution of PopIII stars to the total SFR is very uncertain 
and the parameters of the model will be studied in up- 
c oming work. 

iMcQuinn et al.l ((2007) considered several models for 
the source efficiency, defined to be proportional to M*/M 
for star forming halos of mass M. They considered cases 
where the source efficiency is independent of M or scales 
as Af -2 / 3 or A/ 2 / 3 , but with no redshift dependence. In 
our star formation prescription, the local star formation 



rate is proportional to the gas density and inversely pro- 
portional to the dynamical time (see Eq. [T0|) . When the 
redshift is kept fixed, the source efficiency has negligible 
dependence on mass because the halos have very similar 
density profiles. They are defined to have the same av- 
erage density of 200 times the critical density and have 
similar concentration parameters (see Eq. [J). However, 
for a fixed mass, the source efficiency scales as (1 + z) 3 / 2 
at high redshifts because of the inverse dependence on 
the dynamical time. 

The source efficiency can be a complicated function of 
redshift, mass, and environment. Feedback from pho- 
toionization and supernovae can alt er the SFR locally , 
particularly in the lower mass halos. Illiev et al.l (|2006bf l 
have used a toy model to study the suppression of star 
formation during the reionization epoch. They argue 
that photoheating will raise the Jeans mass scale for grav- 
itational collapse and thereby reduce or halt star forma- 
tion in lower mass halos. Feedback effects are actually 
much more complicated. Photoionization produces elec- 
trons, which are catalysts in H2 formation, and this can 
lead to enhanced cooling. Supernovae will inject energy 
into the ISM and IGM, but the metal enrichment can 
lead to additional cooling. These nonlinear effects are 
best investigated with high re solution , smal l volume hy- 
drodynamic simulations (e.g. iGnedinl 120001 ). They can 
then be incorporated into the star formation prescrip- 
tion for large volume simulations of reionization. 

4.4. Ionization fractions 

We plot the mass and volume weighted ionization frac- 
tions for HI, Hcl, Hell, and electrons in Figure ITD1 In or- 
der to have reionization end at z ~ 6.5, we chose the L50a 
and L50b simulations to have radiation escape fractions 
of 16% and 15%, respectively. Despite the lower escape 
fraction in the latter simulation, there are more ioniza- 
tions at higher redshifts due to the presence of PopIII 
stars. HI and Hel are ionized similarly, both spatially and 
with redshift, because they have similar absorption cross- 
sections and recombination coefficients. Only a small 
fraction ~ 10% — 20% of Hell is ionized by z = 5, with 
more in the L50b simulation because the PopIII stars 
are assumed to have a top-heavy IMF with an effective 
spectrum which peaks at higher frequencies compared to 
PopII stars with a Salpeter IMF. 

In Figure 1111 we show that the computed mass and 
volume weighted residual HI fractions at z < 6.5 are 
both in good agreement with the measurements from 
high redshift quasar s in the SDSS (|Fan et all r2006bh . 
iGnedin fc Fan! (|2006[ ) have shown that high resolution, 
small volume hydrodynamic simulations which resolve 
small-scale objects like Lyman limit systems, can cor- 
rectly count the number of absorptions. Figure 5 in their 
paper shows the agreement in the volume weighted HI 
fraction, but it is not obvious that the mass weighted 
fraction is also in agreement. Their simulations have a 
maximum box length of 8 Mpc/h and will give highly 
biased results for both sources and HII regions. To date, 
no other numerical work involving large volume simula- 
tions have reported this agreement with observations. In 
upcoming work, we will compare the HI optical depths 
(jGunn fc Petersonll965h since it is a more direct observa- 
tional measurement than t he neutral f r action . The mean 
HI fractions calculated by I Fan et aL ( 2006b) are based 




Fig. 10. — Mass (thicker lines) and volume (thinner lines) weighted ionization fractions. In the top panels, we show the HI and Hel 
fractions from simulations with (red, long-dashed) and without (blue, short-dashed) PopIII stars. HI and Hel are ionized similarly, both 
spatially and with redshift, because they have similar absorption and recombination coefficients. In the bottom panels, we show the electron 
and Hell fractions with (orange, long-dashed) and without (green, short-dashed) PopIII stars. Reionization commences earlier and the 
duration is longer with the addition of PopIII stars in the L50b simulation. 



on a model of the density distribution of the IGM by 
iMiralda-Escude et al.l (|2000f ) and this model needs to be 
updated before a precise comparison can be performed. 
Furth ermore, we plan to compar e the dark gap distribu- 
tions ([Paschos fc Normanl [2005) with high redshift ob- 
servations. We can probe the response of the gas to 
reionization by modelling and varying the photoioniza- 
tion feedback in our baryonic prescription. 

4.5. Clumping factors 

In Figure 1121 we compare the clumping factors mea- 
sured at high resolution with the particles and at low 



resolution with the RT grid. The particles can probe 
scales near the PM spacing Axpm = 8.68 kpc/h, but the 
grid is smoothed on scales of Axrt = 278 kpc/h. We 
define the cosmic clumping factors for H and HII as, 



Gh = 



Chii = 



<n H ) 2 



_ \ n HIl) 



<n H > 2 



(21) 



(22) 



where we normalize using the mean cosmic hydrogen den- 
sity (jih)- The clumping factors are calculated using the 
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Fig. 11. — Mass (upper, thicker lines) and volume (lower, thin- 
ner lines) weighted HI fractions with (red, long-dashed) and with- 
out (blue, short-dashed) PopIII stars. The simulated residual HI 
fractions are in good agreeme nt with measurem ents from 19 high 
redshift quasars in the SDSS IIFan et al.H2006bll . 




Fig. 12. — Clumping factors for H and HII measured at high res- 
olution with the particles (thicker) and at low resolution with the 
RT grid (thinner). The njj clumping factors (black) are measured 
using all baryons. The njjn clumping factors are measured using 
ionized gas from simulations with (red, long-dashed) and without 
(blue, short-dashed) PopIII stars. All recombination calculations 
are performed at high resolution using the particles to prevent the 
clumping factors from being underestimated. 



Fig. 13. — The Thomson optical depth from simulations with 
(red, long-dashed) and without (blue, short-dashed) PopIII stars. 
The three horizontal lines are the best-fit (dashed) and 1 — cr 
(dotted) errors from third-year WMA P (r e = 0.088+0.028-0.034; 
IPage et all |2006|; l&pereel et al.ll2007h . The (black, solid) curve 
r e (z) gives the optical depth for instantaneous reionization at red- 
shift 2. 

baryonic prescription which accounts for Jeans smooth- 
ing in identified collapsed halos. Gas clumping factors 
derived using high resolution N-body simulations where 
the baryons are assumed to trace the dark matter at all 
scales are overestimates. 

At high redshifts z > 20 when the density field is still 
highly linear, the H clumping factors from the particles 
and grid are close to unity. As the redshift decreases, the 
abundance of collapsed halos increases and the universe 
becomes clumpier. The values calculated from the parti- 
cles reflect this change, but the grid values are underes- 
timates. At all redshifts prior to complete reionization, 
the HII clumping factors are more severely underesti- 
mated by the grid than for hydrogen. Since HII regions 
originate from highly biased sources with large density 
contrast, the ionized gas will be very clumpy on scales 
much smaller than the grid smoothing scale. By redshift 
z = 6 when reionization is complete, the cosmic clump- 
ing factor is approximately 30, but the grid value is lower 
by a factor of 15. 

Low resolution, large volume simulations have at- 
tempted to account for the subgrid clumping fac- 
tor to improve the recombination c a lculations (e.g. 
Kohler et ail 120051 : llliev et all l2006bl: iMcQuinn etatl 



20071) . For a given density on the RT grid, they have 



applied an averaged correction, but have yet to account 
for the subgrid scatter. We naturally account for scatter 
by directly working with the particles. 

4.6. Optical depth 

The Thomson optical depth to electron scattering up 
to the redshift of recombination z rcc is given as, 



n c (z)-^—dz 

v ; dz 



(23) 



Radiative transfer simulations of cosmic reionization 



13 



where n e (z) is the proper mean electron number den- 
sity at redshift z. WMAP polarization measurements 
with three years of data haye yielded an optical dept h 
r c = 0.088i8;git dPaee et all 120061 : iSpergel et alJl200l . 
which is a factor o f 2 smaller than the first year result 
(|Kogut et al.l l2003). The best-fit value suggests that the 
universe was highly ionized by z ~ 10, but the uncer- 
tainties are still large. 

In Figure [13j we show that the values calculated from 
the L50a and L50b simulations are consistent within 
1 — (7, albeit on the low side. In the latter simulation, 
6 x 10~ 5 of the baryons are turned into PopIII stars by 
redshift z = 6. The addition of this amount of PopIII 
stars has resulted in an increase of Ar c « 0.007 in the 
optical depth. To match the best-fit value would require 
several times more P opIII stars at a fixed total SFR. 
IWvithe fe Ceil (|2006h have argued that the era of PopIII 
star formation can be significantly prolonged if mini- 
halos remain metal-poor until they merge into larger ha- 
los with virial temperatures above 10 4 K where star for- 
mation is assumed to be more efficient. If the redshift of 
complete HII overlap is kept fixed, then a more extended 
reionization epoch can be achieved by introducing more 
PopIII stars at higher redshifts. 

Models of reionization with late overlap at redshift 
Z ~ 6 gene r ally g ive lower values for the optical depth. 
IZahn et"afl (|2007l ) recently obtained r = 0.06 from a 
simulation where reionization ended at z w 6. How- 
ever, their low resolution simulation will systematically 
underestimate the SFR at higher redshifts, resulting in a 
shorter reionization epoch and lower optical d epth. Sim- 
ulation s with earlier overlap as considered bv llliev et al.l 
(|2006al fbh have given larger values which are closer to 
the current best-fit value. However, it has not been 
shown whether simulations with early reionization can 
match observations of quasars at 5 < z < 6.5. Pre- 
viously in §4.41 we discussed that the mass and volume 
weighted residual HI fractions in our simulatio ns are both 
in good agreement with recent observations bv lFan et al.l 
(I2006bft . 

We can account for all halos above the cooling mass 
Af v ir(Tvi r = 10 4 K) for z < 15, but miss relatively more 
of these halos as the cooling mass decreases with increas- 
ing redshift (see Fig. QJ. Consequently, the computed 
optical depths may be underestimated, especially in the 
L50b simulation where PopIII stars make a compara- 
tively larger contribution than PopII stars at higher red- 
shifts. Furthermore, a box size of 50 Mpc/h is not quite 
large enough to calculate the optical depth accurately. 
We have found that small boxes systematically under- 
predict the optical depth. The artificial sudden over- 
lap of the largest HI I region with itself due to periodic 
boundary conditions (Bark aria fc Loeb1l2004D results in a 
less extended reionization epoch and hence a lower value 
for the optical depth. The two test simulations are very 
useful for comparing the differences between reionization 
histories with and without PopIII stars, but larger simu- 
lations are required for detailed quantitative comparisons 
with observations. In ongoing work using a large simu- 
lation with 24 billion particles in a 100 Mpc/h box, we 
find that a modest increase in the SFR at high redshifts 
and in the total PopIII stellar density can give an optical 
depth of r c w 0.09 (Shin et al. 2007). 



4.7. Summary and Conclusions 

We have developed a new hybrid code to run large 
volume, high resolution radiative transfer simulations of 
cosmic reionization. Two test simulations each with 3 
billion particles and 400 million rays in 50 Mpc/h boxes 
have been run to date to study reionization histories with 
and without PopIII stars. Large simulations with 24 bil- 
lion particles 100 Mpc/h boxes will be presented in com- 
panion papers (Shin et al. 2007, Trac et al. 2007). 

The RT calculations utilize an efficient adaptive ray 
tracing algorithm. W e use a ray splitting scheme 
(jAbel fc Wandeltl [2002) to obtain robust angular reso- 
lution for ray tracing from point sources, but introduce 
a new ray merging scheme in order to handle the mil- 
lions of sources. Brute force ray tracing inherently scales 
as 0(N 2 ), but with adaptive merging of co-parallel rays, 
we converge on O(N) scaling as the radiation filling fac- 
tor approaches unity. This ray tracing algorithm can 
be combined with a cosmological hydrodynamic code to 
study the effects of reionization and a nonuniform radi- 
ation field, due to shadowing and shielding, on the high 
redshift Lyman alpha forest. 

A main feature of our RT algorithm, which makes it 
unique compared to others used in large volume reion- 
ization simulations, is that the ionization and recombi- 
nation calculations are performed on the particles rather 
than on a coarse grid. We gain three major advantages 
by working with the particles directly. First, clumping 
factors are calculated at high resolution down to scales of 
several comoving kpc rather than several hundred kpc. 
Second, we account for the time-dependent photoioniza- 
tion cross-sections of neutral gas. Third, we account for 
the self-shielding of radiation sinks. 

We use a particle-multi-mesh N-body code to track 
the gravitational evolution of the collisionless dark mat- 
ter. Collapsed dark matter halos are identified using a 
friends-of-friends algorithm. The high redshift halos are 
strongly biased and linear theory breaks down at smaller 
k than it does for the dark matter density field. The halo 
spectra resemble power-laws, regardless of mass and red- 
shift. We will quantify the halo clustering in more detail 
with the larger simulations. 

The two test simulations can identify halos down to 
a minimum mass of 6 x 10 7 Mq/H. For z < 15, we can 
account for all halos with virial temperatures > 10 4 K 
and some mini-halos, but at higher redshifts, we miss 
some sources and all of the mini-halos. As a result, we 
may be underestimating the reionization process and the 
Thomson optical depth. Un resolved halos can be a dded 
using the method outlined in lMcQuinn et all (|2007h . but 
care must be taken to ensure that input halos at one 
redshift match simulated halos at a lower redshift. They 
have applied the technique to post-processed data, but 
it is more difficult for us since the RT calculations are 
performed concurrently with the N-body calculations. 

Baryons are assumed to trace the dark matter in low 
density regions, but a bias is introduced within collapsed 
halos to account for Jeans smoothing from gas pressure. 
In the present version of the code, the gas is assumed 
to be in hydrostatic equilibrium and follows a polytropic 
equation of state. The derived gas density and temper- 
ature profiles are then used to self-consistently calculate 
the gas clumping factors and star formation rates. We 
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have yet to take into account rotation, cooling, and feed- 
back from photoionization and supernovae, but these ef- 
fects will be incorporated into the baryonic prescription 
in upcoming work. 

The star formation prescription closely resembles those 
used in hydrodynamic simulations. Star formation is re- 
stricted to occur only in halos with virial temperatures 
above 10 4 K and the local SFR is taken to be propor- 
tional to the gas density and inversely proportional to 
the dynamical time. Our cosmic SFR is in good agree- 
ment with hydrodynamic simulations. We find that the 
source efficiency is independent of mass for a fixed red- 
shift, but is redshift dependent when the mass is fixed. 
Photoionization, supernovae, and metal enrichment can 
alter the source efficiency and cosmic SFR. These non- 
linear effects should be investigated with high resolution, 
small volume hydrodynamic simulations and then incor- 
porated into the star formation prescription for large vol- 
ume simulations of reionization. 

We have simulated two reionization models and have 
fixed the redshift of complete HII overlap to occur be- 
tween 6 < z < 6.5. Our results are consistent with 
the latest observations from SDSS and WMAP, both of 
which still have considerable uncertainties. We empha- 
size that our values for the optical depth are for the two 
chosen reionization histories. The actual star formation 
efficiency, number of ionizing photons per baryon, radi- 
ation escape fraction, and the dumpiness of the IGM at 
high rcdshifts are still unknown up to an order of mag- 
nitude. Therefore, it is important to vary the parameter 
space in studies of reionization. This can be efficiently ac- 



Abel, T., Norman, M. L., & Madau, P. 1999, ApJ, 523, 66 
Abel, T., & Wandelt, B. D. 2002, MNRAS, 330, L53 
Barkana, R., & Locb, A. 2001, Phys. Rep., 349, 125 
— . 2004, ApJ, 609, 474 

Becker, R. H., Fan, X., White, R. L., Strauss, M. A., Narayanan, 
V. K., Lupton, R. H., Gunn, J. E., & Annis, J. 2001, AJ, 122, 
2850 

Bryan, G. L., k, Norman, M. L. 1998, ApJ, 495, 80 
Cen, R. 1992, ApJS, 78, 341 
— . 2003a, ApJ, 591, L5 
— . 2003b, ApJ, 591, 12 

Cen, R., & Ostriker, J. P. 1992, ApJ, 399, L113 

Ciardi, B., Ferrara, A., & White, S. D. M. 2003, MNRAS, 344, L7 

Cohn, J. D., & White, M. 2007, ArXiv e-prints, 706 

Dolag, K., Bartclmann, M., Perrotta, F., Baccigalupi, C, 

Moscardini, L., Meneghetti, M., & Tormen, G. 2004, A&A, 416, 

853 

Fan, X., Carilli, C. L., & Keating, B. 2006a, ARA&A, 44, 415 

Fan, X., Narayanan, V. K., Lupton, R. H., Strauss, M. A., Knapp, 
G. R., Becker, R. H., White, R. L._, Pentericci, L., Leggctt, S. K., 
Haiman, Z., Gunn, J. E., Ivezic, Z., Schneider, D. P., Anderson, 
S. F., Brinkmann, J., Bahcall, N. A., Connolly, A. J., Csabai, I., 
Doi, M., Fukugita, M., Geballe, T., Grebel, E. K., Harbeck, D., 
Hennessy, G., Lamb, D. Q., Miknaitis, G., Munn, J. A., Nichol, 
R., Okamura, S., Pier, J. R., Prada, F., Richards, G. T., Szalay, 
A., & York, D. G. 2001, AJ, 122, 2833 

Fan, X., Strauss, M. A., Becker, R. H., White, R. L., Gunn, J. E., 
Knapp, G. R., Richards, G. T., Schneider, D. P., Brinkmann, J., 
& Fukugita, M. 2006b, AJ, 132, 117 

Furlanetto, S. R., Zaldarriaga, M., & Hernquist, L. 2004, ApJ, 613, 
1 

Gao, L., White, S. D. M., Jenkins, A., Frenk, C. S., & Springel, V. 

2005, MNRAS, 363, 379 
Gnedin, N. Y. 2000, ApJ, 535, 530 
Gnedin, N. Y., & Fan, X. 2006, ApJ, 648, 1 



complished by combining numerical and semi-analytical 
methods. 

We will learn more about the epoch of reionization 
from additional observations of high redshift quasars in 
the SDSS and improved measurements of the Thomson 
optical depth from WMAP. The next generation of ob- 
servations, primarily high redshift surveys of Lyman al- 
pha emitting galaxies, CMB measurements of KSZ ef- 
fect from free electrons, and radio observations of the 21 
cm radiation from neutral hydrogen, can provide signifi- 
cantly stronger constraints, but this demands much more 
accurate theoretical calculations. Our approach for large 
volume, high resolution simulations can provide the ac- 
curacy and feasibility required for the task. 
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